/// Script to produce relative risk graph and estimates using nlcom command and coefplot package
/// Fabien Cottier
/// SPEI DATA // quasi_poisson


///////// Graph immediate SPEI effect only(SPEI Y0)

* write programe for relative risk estimation // margins cannot produce rr estimates for cont. variables
program define rr_est	
	local models ""
	local x=-.5
	
	* Loop to create multiple nlcom statement
	forvalues j=1(1)11 {

		if "`2'"=="linear"{
			local models `"`models' (_b[wmean_spei]*`x')"'
		}
		
		else if "`2'"=="quadratic"{
			local models `"`models' (_b[wmean_spei]*`x'+_b[wmean_spei#wmean_spei]*(`x')^2)"'
		}
		
		local x=`x'+1/10
	}

	* estimates relative risk using nlcom
	est resto `1' 
	nlcom `models' , post
	estimates store est_rr
	
	* produce plot using coefplot
	coefplot (est_rr), ///
		eform yscale(range(.5 2.5)) ylabel(.5 "-50%" 1 "0" 1.5 "+50%" 2 "+ 100%" 2.5 "+150%") xlabel(1 "-0.5" 3.5 "-0.25" 6 "0" 8.5 "0.25" 11 "0.5") ///
		ytitle(Relative change in irr. migration (%)) xtitle(Standardized Precipitation Evapotranspiration Index) vertical recast(line) ciopts(recast(rarea)) yline(1) 
end


///////// Split sample graph effects SPEI (SPEI Y0)


* write programe for relative risk estimation // margins cannot produce rr estimates for cont. variables
program define rr_est_Sp	// MuL_Sp for multiple lag with split sample
	local models_L0 ""  
	local x=-.5
	* Loop to create multiple nlcom statement
	forvalues j=1(1)11 {
	
		if "`3'"=="linear"{
			local models_L0 `"`models_L0' (_b[wmean_spei]*`x')"'
		}
	
		else if "`3'"=="quadratic"{
			local models_L0 `"`models_L0' (_b[wmean_spei]*`x'+_b[wmean_spei#wmean_spei]*(`x')^2)"'
		}
		
		local x=`x'+1/10
	}
	
	* estimates relative risk using nlcom // With High labor agriculture
	est resto `1'
	nlcom `models_L0', post
	estimates store est_rr_L0H

		
	* estimates relative risk using nlcom // With Low labor agriculture
	est resto `2'
	nlcom `models_L0', post
	estimates store est_rr_L0L
	
	* graph
	coefplot (est_rr_L0H), bylabel(High Agriculture - Immediate effect) || ///
		(est_rr_L0L), bylabel(Low Agriculture - Immediate effect) ||, ///
		eform yscale(range(.5 2.5)) ylabel(.5 "-50%" 1 "0" 1.5 "+50%" 2 "+ 100%" 2.5 "+150%") xlabel(1 "-0.5" 3.5 "-0.25" 6 "0" 8.5 "0.25" 11 "0.5") ///
		ytitle(Relative change in irr. migration (%)) xtitle(Standardized Precipitation Evapotranspiration Index) vertical recast(line) ///
		ciopts(recast(rarea)) yline(1)  byopts(row(3))
		
end





 
